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ABSTRACT 

We examine the duty cycle and the history of star formation (SFH) for high-redshift 
galaxies at z > 6 using cosmological hydrodynamic simulations. We find that, even 
though individual galaxies have bursty SFH, the averaged SFH between z ~ 15 to 
z = 6 can be characterized well by either an exponentially increasing functional form 
with characteristic time-scales of TOMyr to 200 Myr for galaxies with stellar masses 
Ms ~ lO^M© to > lO^^M© respectively, or by a simple power-law form which exhibits 
a similar mass dependent time-scales. Using the SFH of individual galaxies, we measure 
the duty cycle of star formation (I^Csfh); i-e., the fraction of time a galaxy of a 
particular mass spends above a star formation rate (SFR) threshold which would 
make it observable to the Hubble Space Telescope (HST) during a given epoch. We 
also examine the fraction of galaxies at a given redshift that are brighter than a 
rest-frame UV magnitude (M^v ~ —18), which is sufficient enough to make them 
observable {DCm^^)- We find that both i^CsFH and DCm^^ make a sharp transition 
from zero (for galaxies with Ms ^ 10^ Mq) to unity (for Ms > lO^M©). The measured 
duty cycle is also manifested in the intrinsic scatter in the Ms — SFR relationship 
(~ 1 dex) and Ms — M^v relationship (AM^v ~ ±lmag). We provide analytic fits to 
the DC as a function of Ms using a sigmoid function which can be used to correct for 
catalogue incompleteness. We consider the effects of duty cycle to the observational 
estimate of galaxy stellar mass functions (GSMF) and the star formation rate density 
(SFRD), and find that it results in a much shallower low- mass end slopes of the GSMF 
and a reduction of > 70% of our intrinsic SFRD, making our simulation results more 
compatible with observational estimates. 
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1 INTRODUCTION 

Cosmological studies of structure formation have made 
tremendous progress over the past decade. At low red- 
shift, large galaxy surveys such as the SDSS and 2dF have 
measured th e galaxy distribution with high accuracy (e.g., 
iTegmark eTa l. 2004; Cole et al. 2005). Other cosmological 
probes, such as the distances to th e Type la supernovae 
(iRiess et all Il996l: IPerlmutter et al ] [19991 cluster abun- 
dance ([Bahcall et al. Il997l ), temperat ure anisotropy of cos- 
mic microwave background radiation (Spergel & SteinhardJ 
[2000 : .Komatsu et al...2011. ) are also employed. A combina- 
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tion of these observational results all point to a leading theo- 
retical model of our Universe that is dominated by cold dark 
matter (CDM) and dark energy (A). Our goal is to draw a 
self-consistent picture of galaxy formation and evolution un- 
der the framework of ACDM model. 

The observations of local galaxies often show that the 
SFH of local galaxies can be f itted well with a exp o nentially 
declining!: functional form (Gallaghe r et al.l Il984l : [Sandagq 
ll986l : lKenn icutt 1998). However, in the ACDM model, small 
objects form first, and they later merge into more massive 
systems. Therefore one would expect that there must be 
an early phase of galaxy formation when stars are vigor- 
ously formed in small dark matter halos, and larger galax- 
ies are being assembled for the first time with an increas- 
ing SFH. The so-called Lilly- Mad au diagram ([Lilly et al.l 
Il996l : iMadau et al.l[T99i l reflects this rise and fall of cos- 
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mic SFH. The initial version of this diagram had its peak 
around z ^ 1, however, as the observations became deeper 
and deeper, observers found a large number of fainter galax- 
ies, and soon it was consid ered to be flat d uring z ^ 1 — 3 
or to even higher redshifts (Ste idel et al.|[l99 9') . The redshift 
frontier was pushed to earlier epochs b y the Lyman-break 
drop-out technique (|Steidel et alJIlOOil ), and a larger num- 
ber of galaxies were detected at z = 3 — 6, again showing 
a decline of cosmic star formation rate density (SFRD) to- 
wards highe r-z fe.g.JOuchi et al.ll200i:lBouwens et al.ll20o3 : 
[Ciavalisco et al.ll20oilBouwens et al.ll201lh . 

With the advent of the new WFC3 camera on- 
board the HST, astronomers have begun to catalog large 
num bers of candidate star-forming^alaxies at z > 6 
fe.g. JOuchi" et alj I I2OO9I : iBouwens et al.ll ^id : IWilkin7 et alJ 
I2OI1I : lYan et alJ |201Q| : lOesch et alJ I2OI2I I. This is again 
achieved by the eflicient selection of candidates in the 
color-color plane using the Lyman-break drop-out tech- 
nique. However, without a spectroscopic follow-up, physi- 
cal properties such as stellar mass, SFR, age and metal- 
licity are very difficult to derive accurately. The cur- 
rent analysis relies on the SED fltting of a combination 
of stacked WF C3 rest- frame UV data an d Spitz er IR 
data (e.g.. Labbe et al. '2OIQI: Schaerer de Barro^ |201Q| : 
iMcLure et al.1 I2OIII : iTrenti et al.1 l201lh using population 
synthesis codes. This flt has several free parameters includ- 
ing stellar mass, metallically and age, and it requires an 
assumption about the SFH of the galaxy. Often case, the as- 
sumptions include an instantaneous burst of star formation 
followed by a decline, or a burst with a constant SFR. How- 
ever the SFH is not well understood for galaxies at ^ > 6, 
and it would be helpful to use a physically motivated models 
in deriving the age and stellar mass of these high-z galax- 
ies. As we argued earlier, one would expect an increasing 
mean SFH at high-^, and some researchers have already 
found such SFH in c osmological hydrodynamic simulations 
as well as in some observational data 
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In our previous work, we have measured the evolution 
of galaxy luminosity function (LF) and GSMF in our cosmo- 
logical hydrodynamic simulations (Jaacks et al. 2012). One 
of the notable results was the very steep faint-end (low- 
mass end) slopes of simulated galaxy LF/GSMF, with slopes 
steeper than a=— 2.0at2;>6 and magnitudes fainter 
than Muv = -19 and Ms < 10^ M©. These simulation re- 
sults are somewhat inco mpatible with the current observa - 
tional estimates ( Bouwens et al .11 20 111 : I Gonzalez et al.ll201lh . 
both of whom find less steep faint-ends in the observed LF 
and GSMF. However, at the same time we found that the 
low-mass, star-forming galaxies are necessary to maintain 
the ionization of i n tergalactic medium at z > 6 (see also 
iTrenti et al.1 (I2OIOI): ISalvaterra et al.1 (|201lh : iBouwens et al.1 
(|201lh : I Finkelstein et al .1 (|2012l )). therefore important for 
the reionization of the Universe. Thus we realized that it 
is important to understand the SFH of these high-z galax- 
ies and the duty cycle of star formation in our simulations, 
which is exactly what we will do in the current pa per using 
the same simulation sets as in il jaacks et aP (I2OI2I ). 

The discrepancies mentioned above between our simula- 
tions and recent observation is less evident when looking at 
photometric properties (i.e., color-color diagrams, LF), but 
dramatically more obvious when comparing physical proper- 



ties (e.g., GSMF, SFRD, M,-SFR relation). The difficulty 
of estimating these properties for high-z galaxies is in large 
part due to the intrinsic scatter of derived physical prop- 
erties, which ca n lead to a large scatter in the M/L ratio 
([Bouwens et al .' 2011). By examining the specific details of 
SF in these galaxies, in particular the episodic nature of 
SF history (i.e., duty cycle), we can better understand the 
discrepancies between observation and theory. 

The idea of using the SF duty cycle as an explana- 
tion for the offset between the observed and theoretical 
GSMF has been demonstrated in previous works. For ex- 
ample, I Lee et ahl (| 20 111 ) considered a scenario in which the 
SFHs of lower luminosity galaxies are bursty. They used the 
observed Ms — Muv relation at 2; = 4 and it's deviation from 
the smooth growth model to estimate the duty cycle. The 
model with DC successfully lowers the number density of 
low-mass objects that would be observed. The goal of this 
paper is to provide the measurement of the duty cycle of 
star formation of simulated galaxies in our cosmological hy- 
drodynamic simulation without any ad hoc assumptions on 
gas dynamics. 

This paper is organized as follows. We briefiy describe 
our simulations in Section[2l We then present our results of 
SFH and duty cycle in Section [3Jl 13.21 We discuss the im- 
plications of our results on our simulated GSMF and SFRD 
in Section[4l and summarize and discuss in SectionjS] 



2 SIMULATIONS 

We use a modified version of the smoothed particle hy- 
dro d^m^micsJ^PH) code GADGET-3 (originally described 
m ISpringell l2005l ). This modified code includes radia- 
tive cooling by H, He, and metals (jChoi Nagamin3 
I2OO9I ) , h eating by a un if orm UV backgr ound of a 
modified Haardt & Madau ("1996') spect rum (jKatz et al.l 
1996 : Dave et al. 1999 : Faucher-Giguere et al.1 l2009l ). an 



lEisenstein Hul (|l999l ) initial po wer spectrum, star for- 
matio n via the "Pressure mod el" (jSchave &: Dalla Vecchial 
2008; IChoi k Nagamind I2OIOI ). supe rnova feedback, sub- 
resolution model of multiphase ISM (Spr ingel HernquistI 
2003 ). and t he multicomponent var iable velocity (MVV) 
wind model (jChoi Nagamind 1 20 111 ). Our current simula- 
tions do not include AGN feedback. 

We have shown earlier that the metal line cooling en- 
hances star formation acr oss all redshifts by about 10 — 30% 
(|Choi Nagamind [2OO9I ). and that the Pressure SF model 
suppresses star formation at high-redshift due to a higher 
threshold density for SF fChoi & Nagamine" 2010) with re- 
spect to the earlier model by Springel & Hernquist (200^). 
IChoi k Nagamind (|201ll ) also showed that the MVV wind 
model, which is based on the momentum-driven wind, 
makes the faint-end slope of GSMF slightly shallower com- 
pared to the constant ve locity galactic wind model of 
Spri ngel k Hernauisti (|2003l ). 

Simulations are setup with either 2 x 400^ or 2 x 600^ 
particles for both gas and dark matter. Multiple runs were 
made with comoving box sizes of 10/i~^Mpc, 34/i~^Mpc 
and 100/i~^Mpc to cover a wide range of halo and galaxy 
masses. Throughout this work they will be referred to as 
N400L10, N400L34 and N600L100 runs, respectively. Sim- 
ulation parameters are summarized in Table [U and the 
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Run 


Box size 


Np 






e 






(DM, Gas) 


(/i-iMo) 


(/z-iMo) 


(/i-^kpc) 


N400L10 


10.0 


2x400^ 


9.37x10^ 


1.91x10^ 


1.0 


N400L34 


33.75 


2x400^ 


3.60x10'^ 


7.34x10^ 


3.38 


N600L100 


100.0 


2x600^ 


2.78x10^ 


5.65x10'^ 


4.30 



Table 1. Simulation parameters used in this paper. The parameter Np is the number of gas and dark matter particles; moM and 
are the particle masses of dark matter and gas; e is the comoving gravitational softening length. 



adopted cosmological parameters are based on the WMAP 
data: 0.26, VLk 0.74, .044, h 0.72, 

(78 = 0.80, Us = 0.96 (iKomatsu et al.ll2009l l. 

At each time step of the simulation, the gas particles 
that exceed the SF threshold density (nfjf = 0.6 cm~^) 
are allowed to spawn star particles with the probability 
consistent with the intended SFR. Each star particle is 
ta^^ed by stellar mass , metallicity, and formation time. 
IChoi Nagamind (fioiol ) have shown that the Pressure SF 
model with the above n^if produces fav orable results com- 
pared to the observed iKennicuttI (| 19981 ) law, in particular 
at the low column density end. Galaxies are then identified 
by a group-finder algor ithm usin^ a simplifi ed variant of the 
SUBFIND algorithm (|SDringel et al.ll200lh . Collections of 
star and gas particles are grouped based on the baryonic 
density field. Properties such as stellar mass, SFR, metallic- 
ity and position, among others are recorded in galaxy prop- 
erty files. A more deta iled description of the group-finder 
process can be found in lNagamine et aP (| 20041 ). 

2.1 Numerical resolution 

To demonstrate that the simulated Ms — SFR relationship 
which is used extensively in the work is resolution indepen- 
dent, we present a convergence study in Fig. [1] Runs with 
2 X 144^ (magenta circles), 2 x 216^ (cyan triangles) and 
2 x 400^ (grey squares) all in a comoving 10/i~^ Mpc box 
size are compared to the full dynamic range median relation- 
ship (black solid line, see Fig. [4|) at z = 6. The results from 
the three runs with different resolution is consistent with 
the black solid line, suggesting that our result on Ms — SFR 
relationship is not affected by the resolution. 

The vertical dashed black line represents our imposed 
resolution limit for the N400L10 run and corresponds to a 
galaxy with approximately 66 star particles. Galaxies below 
this resolution limit are excluded from consideration when 
calculating the Ms - SFR, Ms - M^v, GSMF, k SFRD. 
I Jaacks et aP (|2012) have shown that the gas particle mass 
of the N400L10 run is lower than the typical Jeans mass at 
z = 6 by a factor of ^ 1 — 100, but the run does not satisfy the 
iBate BurkertI (| 19971 ) mass resolution criteria. Therefore 
our simulations are not resolving the collapse of star-forming 
molecular clouds directly, which is why we need to employ 
the sub-resolution multiphase model for star formation in 
our current simulations. 



3 RESULTS 

3.1 Star Formation History 

Fig. [2] shows the mean SFH for all galaxies found in our 
three runs, N400L10, N400L34 and N600L100. It is com- 




logM, [MJ 



Figure 1. Ms — SFR relationship at z = 6 from different reso- 
lution run, 2 X 144^ (magenta circles), 2 x 216^ (cyan triangles) 
and 2 x 400^ (grey squares) all in a 10/i~^ Mpc box size. The 
solid black line represent the full dynamic range (Fig. |4]) median 
fit to this relationship which go as Ms oc SFR^-^^. The vertical 
dashed black line represents our imposed resolution limit for the 
N400L10. 

puted as follows. First, the SFH for each identified galaxy 
is calculated by counting the number of star particles that 
were created in each time bin of At = 10 Myr. Galaxies are 
then grouped by their stellar masses at z = 6, and the mean 
SFH is plotted for each stellar mass group. Note that the 
SFH at higher redshifts includes all the star formation in 
all progenitor galaxies, which may have had lower stellar 
masses than the final mass at z — 6. 

To quantify the mean SFH, we first perform a least 
square fit to an exponential function (solid red line in Fig. [2]), 

SFR = A exp [Me Gyr-'], (1) 

where t is the cosmic time since the Big Bang, A is the 
normalization at z = 6 in units of MoGyr"^, tz6 is the 
cosmic time at 2; = 6, and te is the characteristic time-scale. 
We vary the values of A and te to obtain the best-fitting 
function. We find that the mean SFH can be described fairly 
well by the increasing exponential function for all stellar 
mass groups with a te value ranging from 70 to 200Myrs 
(see Table[2]for the best-fit parameters). This procedure was 
repeated for galaxies at z = 7, 8, 9 with similar results (not 
shown here). 

We also perform the same least square fitting procedure 
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IVLciss RcingG 


A 








(MoGyr-i) 


(Gyr) 




Ms > W^Mq 


76.70 


0.21 


0.9212 


10^ < Ms/Mq < IQio 


12.30 


0.16 


0.9700 


10^ < Ms/Mq < 10^ 


3.17 


0.11 


0.9976 


10'^ < Ms/Mq < 10^ 


0.42 


0.10 


0.9990 



Table 2. Best-fitting parameter values to the exponential func- 
tion (Eq.[T} for galaxies in different stellar mass ranges, obtained 
by the least squares fit of the star formation histories at z > 6. 
The value of A is the normalization constant, te is the character- 
istic exponential time-scale, and is the correlation coefficient. 




Time [Gyr] Time [Gyr] 



Mass Range {SFR) a 

{h-^MQ) (M0 Gyr-i) (Gyr) 

M^ > 10^^ Mq 15.24 059 3^05 0.9592 

10^ < Ms/Mq < 10^0 1.88 0.62 3.74 0.9886 

10^ < Ms/Mq < 10^ 0.30 0.66 5.12 0.9949 

10'^ < Ms/Mq < 10^ 0^05 0.73 8.55 0.9748 

Table 3. Same as Tabled but for the power-law function (Eq.[2]). 
The SFR takes the value of (SFR) at t = rp, and a is the power- 
index of the power-law function. 



Figure 2. Average star formation histories for different groups 
of galaxies divided by their stellar mass at z = 6 with time bins 
of At = lOMyr. The red solid lines show the least square fit to an 
exponential function, and the blue dashed lines show the power- 
law fit. The magenta dashed-d otted line represents th e polynomial 
fit to SFR/ {SFR) found bv iFinlator et al.l (l201lh . To make a 
valid comparison the value of (SFR) is computed for each panel 
and used to show the function. The total number (N) of galaxies 
found in each mass group is also shown in the plot. The best-fit 
parameters can be found in Tables [2] & [3] 



to a power- law function (blue dashed), 



SFR 



= [MoGyr-^], 



(SFR) ~ \tp) '^'^^ ^' 

where (SFR) is the mean SFR calculated between t = 
and t = t^e, tp is the characteristic time-scale when SFR = 
(SFR), and a is the power index. See Table [3] for the best- 
fitting values. 

In the last column of Tables [2] and O the values of corre- 
lation coefficients are listed as a measure of goodness of fit. 
A value closer to unity indicates a better fit. In most cases 
the power-law function does a slightly better job of fitting 
the shape of the SFH, but they both provide a very good 
fit. 

iHernquist Springell (|2003l l found that the cosmic SFH 
in their cosmological SFH simulations can be fitted well by 
oc exp(— 2;/3) at high-redshifts, and speculated that this 
exponential behavior is related to the growth of dark matter 
halo mass function. Our results support their argument, and 
we also suggest that the early star formation history is driven 
by the gravitational instability of dark matter structure, and 
it is natural to expect an exponential or increasing power- 
law behavior of SFH as a function of cosmic time. 

Included in Fig. [2] is a co mparison to a polyno mial fit 
of the average SFH found bv IFinlator et alJ (|201lh . which 
is indicated by the magenta dashed-dotted line. Using the 
GADGET-2 SFH code, they also found increasing SFH, and 



argued that the shape of SFH to be mass invariant and only 
changes by a scale factor through different mass bins, with 
the following functional form: 



3.9r 



(3) 



for t > 2.7 and equal to zero elsewhere, where t is the age 
of the Universe in Gyr. One can see that our simulation 
results do not agree well with this functional form. In con- 
trast, in our simulations we find that both the exponential 
and power-law fits to have an evolving r with varying galaxy 
stellar mass. The source of this difference might be the differ- 
ent models of star formation and feedback in the two codes. 

In addition to the mean SFH, we also examined the 
SFH of individual galaxies as shown in Fig. [3] It shows 
that, although the mean SFH may be smoothly rising, the 
SFH of the individual galaxies are dominated by intermit- 
tent bursts throughout their evolution up to 2; = 6. These 



burst are thoug ht to be the result of a cornbination of galaxy 
merger events (iMihos k Hernauisti ITqQgI : iKennicuttiTToQi : 



K 

iKhochfar Sillj I2OII) and cold mode accretion onto the 
galaxies ("Kere s et al.l 2005). Note that the SFHs shown in 
Fig. [3] include all the star formation in all progenitors that 
end up in the galaxy at z = 6. We now quantify the bursty 
nature of star formation in terms of duty cycle in the next 
section. 



3.2 Duty Cycle 

The definition of duty cycle (DC) of star formation can be 
somewhat ambiguous. In this work we define it in two differ- 
ent ways. The first, DCsfu is based on the SFH of a galaxy, 
and is characterized by the fractional time that a galaxy 
spends above a threshold SFR which is sufl^icient to make it 
observable. The second, DCm^^ examines the number frac- 
tion of galaxies (in a mass bin at a given instant) that are 
brighter than the threshold rest-frame UV magnitude M^v, 
making it observable. 
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Figure 3. Star formation histories of individual galaxies with 
different stellar masses as indicated in each panel taken from the 
N400L34 run, with time bins of At = 10 Myr. On average, they 
all increase gradually, but individual galaxies show intermittent 
bursts of star formation. 



The SFR threshold (SFRth) used to determine DCsfh 
is obtained by converting the HST's magnitude limit of 



MAB,lin 

mula, 



29 to the rest-frame M^v via the standard for- 



Mab,] 



5 log 



lOpc 



+ 2.51og(l + z). (4) 



At z = 6,7,8, it gives M^J ^ -17.74,-17.98,-18.19, re- 
spectively. Then, M^v can be converted to SFRth via the 
relationship between Muv and SFR of galaxies found in our 
simulations which goes on average as SFR oc M^°"^° (Fig- 
Ill). This results in SFRth = 0.50 M© yr~-^ with no consider- 
ation for dust extinction, and SFRth = 1.20 Moyr"^ when 
including Eb-v — 0.10. This value of Eb-v is chosen to be 
consistent with the value used to match the observed rest- 
frame UV luminos ity function (LF) in our previous work 
(|jaacks et al.ll2012l ). and it is centered between the follow- 
ing two recent observations: iBouwens et al.l (|201lh argued 
for little to no extinction at the f aint-end of the LF at 
z = 6, whereas IWillott et ai] (|2012l ) found a best-fit value 
of Av = 0.75, which corresponds to Eb-v ^ 0.19 assuming 
Ry 4.05 (Calzetti et al. 2000l) at the bright-end of the 
UV LF at ^ = 6. This moderate a mount of extinction is 
also c o nsistent with the estima t es by Schaerer de Barroj 
(l201Ql ): EI Barros et al ] (|2012l ): ISchaerer et al ] (I2OI2I ) who 
included nebular emission lines in their SED fits. Since the 
conversion between and SFRth depends on the rela- 

tionship between Muv and SFR in the simulation, the value 
of SFRth depends on the assumed value of Eb-v- 

In Fig. [3] we plot unc orrected observationa l estimates 
based on SED fits found in de Barros et al.l (|2012l V Although 
blue triangles and squares are median points from the same 
data set they are derived through SED fits which make dif- 
ferent assumptions regarding SFH. The blue triangles are 



calculated using a rising SFH, whereas the squares are cal- 
culated using a constant SFH. Both include nebular emis- 
sion lines and Lya. At first look it would seem that our 
data (yellow diamonds, squares and triangles) have better 
agreement with the blue squares which assume a constant 
SFH. This however is slig htly misleading due to the fact 
that (|de Barros et al.|[2012 ) find that when using the rising 
SFH on average a high extinction is also obtained. Therefore 
this comparison is completely dependent upon our included 
amount of extinction, i.e. if we were to include more extinc- 
tion we would have better agreement with the blue triangles. 
This is further justification to explore a range of extinction 
as we do in this work. 

Recent 2; > 6 publications teouwens et al. I I2OIII : 
[Smit et al.ll2012l ) suggest a Muv-dependent extinction based 
on an observed relationship between Muv and the UV- 
continuum slope (/3). While these results are interesting, we 
have chosen not to adopt this model, and instead apply con- 
stant extinction in an effort to show a range of Eb-v values 
and their impact systematically. Any results obtained from 
using a Muv-dependent model of extinction with fall within 
the range presented in this work. 

We find that the size of time bins Atduty for computing 
SFH makes a difference in determining the value of DC. We 
determine its value based on the mass resolution of each 
simulation and SFRth, i.e., Atduty = nistar/ SFRth, where 
mstar is the mass of star particles in the simulation, which 
is taken to be the half of gas particle mass in the simulation 
setup. This method insures that the simulation is able to 
satisfy the threshold SFR if one star particle forms in any 
given time bin and allow for the comparison between runs 
with different mass resolutions. 

Based on the above definition, we adopt time bins 
of Atduty = 10 Myr and 0.27 Myr for the N400L34 and 
N400L10 runs, respectively, for no extinction, and Atduty = 
4 Myr and 0.22 Myr for the case with Eb-v = 0.10. The 
system time steps of these simulations are typically Atsys = 
10^ — 10^ yrs at z = 6 — 7, therefore we have at least 10 
system time steps in each Atduty, providing reasonable time 
resolution in each Atduty 

Three example SFHs between z = 6 — 7 are shown in 
Fig. [5] for galaxies in the N400L34 run with stellar masses 
ranging ^ 10^ - lO^^M© and DCsfu = 0.08, 0.48, 1.00. This 
figure illustrates the method used to calculate DCsfh, which 
is the fraction of time spent above SFRth = 1.20 Moyr"^ 
(red dashed line). 



3.2.1 Duty cycle vs. Galaxy Stellar Mass 

As we described above, we compute both DCsfh and 
DCm^^ for each galaxy in each simulation. We then make 
a scatter plot of duty cycle as a function of galaxy stellar 
mass Ms, bin the data in each log Ms bin, and obtain the 
median in each bin. Figure [6^, b,c show the results of DCsfh 
at 2; = 6, 7, 8 with dust (red triangles), without dust (blue 
circles), and DCmuv with dust (dashed black line), without 
dust (solid black line). The top axes of panels (a), {h) and (c) 
show the number density of objects (t){Ms) at a given mass 
taken from the Schech ter fits to the GSMF at each redshift 
found in Jaacks et al.l (|2012). 

We find that both methods of DC exhibit a characteris- 
tic steep transition from DC = to unity, and pass through 
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Figure 6. Duty cycles {DCsfh and DCm^^) for galaxies in the N400L34 and N400L10 runs at z = 6, 7, 8. Each method is calculated with 
(red triangles, dashed line) and without dust extinction (blue circles, solid line). The bottom right panel shows the redshift evolution of 
^C'sFH without dust extinction. The dotted line has been added to indicate DC = 0.50 to aid comparison, and error bars represent one 
standard deviation in each mass bin. The top axes of panels (a), (6) and (c) show the logarithm of comoving number density {\og(f){Ms)) 
of objects at a given mass taken from Schechter fits to the intrinsic GSMF at each redshift found in Jaack s et al.l (|2012l ). The blue and 
red lines going through the data points represent the fit to Equation [5] for each of the data sets. Fitting parameter values can be found 
in Table [5] 



DC = 0.50 at approximately the same mass bin at each 
redshift. See Table g] for the mass-scales at which the DC 
crosses the value of 0.50. This transition can be modeled 
well by employing a sigmoid function described by 



DC{Ms) = 



exp 



log(M, 



+ 1 



(5) 



where due to its form the parameter a corresponds to the 
value of log (Ms) where DC = 0.50, and 6 is a measure of 
the steepness of the transition; i.e as b ^ the function 
approaches a step function, and as 6 ^ oo the transition 
becomes infinitely flat. 

Least square fits to Equation [5] were performed for 
each DC at each redshift and the results are shown in 
Fig. [6^,b,c,d as the blue, red, black and black dashed lines 



in each panel. Values of the fit parameters for each model 
at each redshift can be found in Table \5\ The values of a 
parameter is essentially the same as those given in Table [H 
and the closeness of the a values in the two tables show that 
the sigmoid function provides a good fit. 

We note that the N600L100 run result was omitted from 
Fig. [6l because its value of Atduty (=30 Myr) exceeded the 
dynamical time of typical systems at these redshifts, which is 
about 10 Myr for a galaxy of Ms = 10^ Mq. In other words, 
the N600L100 run has insufficient mass resolution to esti- 
mate the duty cycle at these redshifts reliably. 

As shown in Fig. [6]i, we find that the transition mass- 
scale shifts to a lower mass at higher redshift, simply re- 
flecting the fact that the galaxies at higher redshifts are less 
massive than at lower redshifts. Similar evolution is found 
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Figure 4. Relationship between Muv — SFR at z = 6 with 
Eb-v — 0.10. Yellow squares, triangles and diamonds repre- 
sent N400L34, N600L10 and N400L10, respectively. The black 
open triangles are the median found in each bin with error bars 
representing one standard deviation. The solid black line is the 
least square fit to the median data poi nts. The blue tr i angles 
and squares represent observations from Ide Barros et al.l (|2012l ) 
where the triangles are based on SED fits which include rising 
SFH with nebular emission lines and Lya and the squares are 
based on SEDs with constant SFH with nebular emission lines 
and Lya. 
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Figure 5. Examples of three different SFH computed with 
^^duty = 4 Myr from our N400L34 run between z = 6 - 7 and 
corresponding i^CsFH for galaxies of varying masses. The thresh- 
old of SFRth = 1.20 M0 yr"-*^ used to calculate DCsfh is shown 
as the red dashed line. 



in DCmuv ? which is not shown here. This redshift evolution 
trend might seem counter-intuitive at first, because it sug- 
gests that we would see more lower mass objects at 2; = 8 
with high DCs than at 2; = 6. However by examining the 
DC as a function of number density (Fig.[6l top axes) we can 
clearly see that, for a given (j){Ms), the DC at z = 6 is higher 
than at 2; = 8. For example, for galaxies with a comoving 
number density of (j){Ms) — 10~^ [Mpc~^], our results give 
i^CsFH ^ 0.20, 0.15, 0.08 at z = 6, 7, 8, respectively. 

Note th at our I^C s fh d efinition is similar to the one 
put forth bv iLee et all ([2Q0a ). who defined the DC to be 
the typical duration of star formation with respect to time 
span covered by the survey. Since their work is at lower 
redshifts than ours, making direct comparisons is not really 
appropriate. However, the two results seem to be roughly 
consistent with each other as they find that galaxies with 
Ms < 10^ '^ M© are mostly UV-faint and have DC = 0.2-0.4 
at ^ ^ 4 - 5. 

Although individual galaxies in our simulations do ex- 
hibit bursty SFHs, this is by no means the only scenario 
in which a galaxy can exhibit a DC. It is also possible for 
galaxies with rising or declining SFHs to have a DC. Take 
for example the average SFHs broken down by stellar mass 
presented in Fig.[2l If we use Equation ([1]), used to describe 
our exponential rising mean SFH, and apply the concept of 
DCsFU during the time period between z = 6 — 7 we find 
that the DCsfu ~ 1.0,1.0,0.61,0.0 for the case with dust 
and DCsFH ~ 1.0,1.0,1.0,0.0 for the case without. These 
transitions from DCsfu = 1.0 to DCsffl = 0.0 occur in the 
same mass range as those calculated above. 

This suggests that a certain DC value does not corre- 
spond to a unique SFH, and the mapping between DC and 
SFH is certainly not one-to-one. However we consider that 
it is unlikely that all galaxies have a smoothly rising SFH 
with the same functional form, therefore the bursty SFH 
would be more natural outcome of a hierarchical structure 
formation under the ACDM paradigm. 



3.2.2 Future Observations 

Since both DCsfu and DCm^^ are highly sensitive to the 
AB magnitude limit of the observing instrument it is an 
interesting exercise to look a the effect of increase sensitiv- 
ity. The James Web Space Telescope { J WST) is planne d to 
have a magnitude limit of MAB,iim ~ 31 ((Gardner et al.l 
I2OO6I ) which, using equation Q, would equate to a M^v ~ 
— 15.75 at z = 6. With this this new threshold, we cal- 
culate a SFRth = 0.08 [M© yr"^] with Eb-v = 0.0 and 
SFRth = 0.30 [Mq yr"^] with Eb-v = 0.10, both of which 
are well below the values calculated above. Lowering both 
the M^v and SFRth leads to a shift in the DC transition to- 
wards lower masses, and this new transition mass-scales can 
be seen by the dashed black lines in Fig. [71 For the JWST 
magnitude limit, the transition galaxy stellar mass-scale is 
lower by about 0.5 dex and 0.8 dex for Eb-v = 0.0 and 
0.10, respectively. The JWST should be able to see galaxies 
down to Ms ^ WMq at z = 8. 
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Figure 7. The transition stellar mass-scale where DCm^v ~ ^-^^ 
for both Eb-v = 0-00 (blue circles) and Eb-v = 0-10 (red trian- 
gles) calculated using HSVs AB magnitude limit (Mj^B,lim ~ 29) 
shown by the solid black line and JWSVs limit {MAB,lim ~ 31, 
dashed black line). 
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0.00 


7.61 




DCsFH 


0.00 
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Table 4. The values of log(Ms/M0) at which DCm^^ and 
DC^yh cross 0.50 at 2; = 6 — 8. 
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0.00 
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0.00 


7.62 


0.158 




DCsFH 


0.00 


7.63 


0.309 






0.10 


8.08 


0.136 




E>CsFH 


0.10 


8.10 


0.331 



Table 5. Best fit parameter values of the sigmoid function 
(Eq. [5]) used to characterize DC. The Eb-v column shows the 
assumed value of extinction in determining the M^ and S'Fi^th 
(see Eq. |4}. 



4 APPLICATION OF DUTY CYCLE 

Fundamentally DCm^^ and DCsfh are measure of the in- 
trinsic scatter in the relationships between Ms and Muv 
(Fig. [8]), or between Ms and SFR. For each of the rela- 
tionships, we find the scatter of about AMuv ±1 mag 
and ~ 1 dex in SFR, respectively (Fig. [9]). The duty cy- 
cle must be considered when estimating the intrinsic values 
from the observed results, and when making comparisons 
between simulations and observations. 

In Fig. [8] we show the Ms — Muv relationship at z = 6 
with Eb-v = 0.10 for N400L10 (yellow diamonds), N400L34 
(yellow squares) and N600L100 (yellow trian gles) runs, plot- 
ted wi th median points f rom observations b y lGonzalez et alj 
(| 20 111 blue squares) and lStark et alJ (|2009l , red triangles) as 
well as 13 spectroscopi cally confirmed z — 5.5 — 6.5 galax- 
ies presented in I Curtis- Lake et al.l (|20 12, black cross es). The 
dashed blue line is a ^ = 4 fit from the same Gonzalez et alJ 
(|201lh work as above. Here we find good agreement be- 
tween observations and our simulations, although our me- 
dian trend (solid black line) deviates from the z — A Gon- 
zalez relationship toward the low-mass end. Further obser- 
vations will be needed at lower masses to determine if this 
deviation is significant. 

In Fig. [9] we demonstrate the scatter in the Ms — 
SFR relationship with each o f our runs a t z — 7, 
shown with observat i ons from iLabbe et al ] (|2010h and 
ISchaerer k, de BarrosI (|2010l V Agreement with observations 
is somewhat tentative as the scatter in the observed points 
is larger than in the simulated data. This, in part, can be at- 
tributed to the debate as to whether or not nebular emission 
lines ( NELs) need to be included in t he SED fittings at high- 
z. The ISchaerer &: de BarrosI (|2010l ) da ta points (blue) are 
obtained with NELs included, while the lLabbe et al.l (|2010h 
points (red) do not include NELs. 

Currently there are inconsistencies between observa- 
tions and theoretical predictions of the GSMF and SFRD. 
Simulations predict steeper low-mass slope s of the GSMF 
(om) than current observational estimates ([Gonzalez et alJ 
l201lh with a difference of Aq^m ^ 0.80 at z = 6. This dif- 
ference is more significant than that of the faint-end slope 
of the UV LFs, wh ich is A^l ^ 0.40 (|Bouwens et al.ll201ll : 
I Jaacks et alJlioi^ ). Likewise the theoretical predictions of 
SFRD at z = 6 — 10 are in conflict with the estimates ob- 
tained by observations, showing a n offset of more than one 
order of magnitude in total SFRD. IChoi Nagamind (|2012l ) 
showed that these discrepancies can be partly understood if 
we consider the mass and flux limits of the current high- 
z galaxy surveys, however, they did not consider the effect 
of the duty cycle. In the following sections we show that, 
through the application of the DC ^ we can converge obser- 
vations and theory to make a more consistent portrayal of 
the early Universe. 



4.1 Galaxy Stellar Mass Function 

The observed GSMF is estimated using a combination of 
the observed UV LF and Muv — Ms relation, both of which 
are quantities that have the DC built in (i.e., what can- 
not be observed cannot be included). Our simulations give 
the intrinsic number of galaxies in each mass bin, and this 
does not take into account of the observable threshold in 
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Figure 10. GSMF (solid lines) at 2; = 6, 7, 8 with both DCsfh (dashed lines) and DCm^^ (dot-dashed line) applied. This procedure is 
done with the effects of dust extinction (ri ght column ) and without (left column) at each redshift. Yellow shade represents the GSMF 
derived from observations presented in Gonzal ez et al.l (|201lh . The crosses in panels a,b represent simple Muv cut of Muv = — 18 applied 
to the each run prior to computing the GSMF. 



the rest-frame UV light. Therefore we have to apply the ob- 
servational threshold captured by the DC before a proper 
comparison with observations can be made. The application 
of DC as a function of Ms reduces the number of galaxies 
in certain mass bins where DC < 1. 

In Fig. [To] we explore the effect of DCsfh (dashed line) 



and DCm^^ (dot-dashed line) on the GSMF. The solid lines 
in all panels represent the GSMF found in ou r simulations 
and originally presented in Jaack s et al. 1 (2012). The yellow 
sh aded areas repre s ent th e GSMF derived from observations 
bv lGonzalez et al.1 (|201lh . who utilize SED fits to obtain the 
Ms — Muv relationship and then a combination of the UV LF 
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Figure 8. Relationship between Ms — Muv at z = 6. Yellow 
squares, triangles and diamonds represent data from N400L34, 
N600L100 and N400L10, respectively. Black open squares are 
the median found in each bin with the error bars repre- 
senting one standard deviation. Solid black line is the least 
square fit to the median data points. Blue circles and red 
triangles re present median points of uncorrected ob servations 
taken from loonzalez et al.l (|201ll ): IStark et all (|2009l ) and the 
black crosses are f rom r ecent CANDELS observations found in 
ICurtis-Lake et al.l (|2012l ) . The blue d ashed line is a fit to 2; = 4 
data points also found in Gonzalez et al. (2011); note that this is 
not a fit to the z = 6 data set presented. A horizontal dashed black 
line represents the approximate UV magnitude limit of HST. 



and the mass-to-light to determine the GSMF. We present 
these results with (right column) and without (left column) 
extinction correction as we described in Section 13.21 These 
results are obtained through 

0obs(^, Ms) = 0int(^, Ms) X DC{z, Ms), (6) 

where (/)obs(^, Ms) is the observed number density at a given 
redshift and stellar mass Ms, and (/)int(^,Ms) is the value 
intrinsic to our simulations at the same redshift. 

Taking the DC effect into account results in the slope 
of the low-mass end of GSMF becoming shallower at all red- 
shifts by AttM ~ 0.26, where Aqm — |aM(obs) — aM(int)|. 
This decrease of aM brings our GSMF into better agree- 
ment with observations at the low-mass end. The difference 
seen between the left column {Eb-v = 0.0) and right col- 
umn (Eb-v = 0.10) in Fig. 1101 can be attributed to the fact 
that the DCs calculated with dust extinction tend to tran- 
sition from zero to unity in a mass range which is ^ 0.50 
dex higher than those without. This is an expected effect, 
because as extinction increases less star formation will be 
visible, resulting in a lower duty cycle at a given mass. 

In terms of faint-end slope reduction, the application of 
both DCmuv sind DCsfu obtain similar results to the use of 
simple flux limit applied to our galaxy population, as seen 
by the black crosses in Fig. [TOl Deviations between DCm^^ 
and a simple flux limit can be attributed to the fact that the 
the DC is a fit to a scatter in a particular mass bin (error 



Figure 9. Relationship between Ms — SFR at z = 7. Yellow 
squares, triangles and diamonds represent data from N400L34, 
N600L100 and N400L10, respectively. Black open squares are the 
median found in each bin with the error bars representing one 
standard deviation. Solid black line is the least square fit to the 
median data points. This fit is also representative of the average 
relationship at z = 6,7,8, where SFR oc M^-^^. Blue squares 
and red circles show HST/IRAC detected objects at z = 7 taken 
from Schae rer de Barrosi (|2010); Labb e et al.l (|2010|), while the 
blue tria ngles represent non IRAC detected objects at z = 7 
(|Schaerer de Barro^l2010l ). 



bars in Fig. [6]), thereby allowing for the inclusion of more 
objects which would otherwise be excluded by a simple cut. 

4.2 Star Formation Rate Density (SFRD) 

The SFRD is also highly effected by the DC . The observed 
SFRD is obtained by integrating the UV LF to the observed 
flux limits (Muv — 18 at ^ = 6) and converted to SFRD 
via the standard conversion factor depending on the IMF. 
Therefore if we wish to compare our simulated SFRD to the 
observed SFRD, we must take the DC into account. For 
example, if only a fraction of galaxies of a given Ms are ob- 
servable, then only that fraction can be allowed to contribute 
to the total SFRD at a particular redshift. In principle this 
is similar to applying a flux limit to our simulated galaxies 
and the exclusion of those which do not reach an observable 
limit, and the same effect can be taken care of by the DC . 

Fig. [11] demonstrates the effect of applying the DCsfh 
with Eb-v = 0.10 to the SFRD at z = 6, 7, 8. By breaking 
down the galaxy population at each redshift into mass bins 
of 10^-^ < Ms/Mq < 10^ (dashed red), 10^ < Ms/Mq < 10^ 
(dot-dashed blue) and Ms > lO^M© (dotted black), we are 
able to examine the contribution of each to the total SFRD. 
This also allows us to apply our DCsfu to each of galax- 
ies contained in these mass bins and remove those galaxies 
which are not currently observable from our total simulated 
galaxy population. To account for the rapid transition in the 
DCsFU and the large bin size of the decomposed SFRD, we 
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Figure 11. Star formation rate density (SFRD), both intrin- 
sic (SFRDint) and derived observable (SFRDobs)? as a function 
of redshift shown with DCsfh apphed. The result for different 
galaxy stellar mass ranges are shown: 10^'^ < Ms/Mq < 10^ 
(dashed red), 10^ < Ms/Mq < 10^ (dash-dot blue) and Ms > 
10^ Mq (dotted black). The total observable amount of all mass 
ranges (SFRDofes) is shown as the dark green solid line. Total 
SFRDint without DCsfh applied (i.e. the true total in our simu- 
lation) is shown by solid magenta line with crosses. Observational 
estimates taken fro mfgouwens et al. (2011) (black squares). See 
Ijaacks et al ] (l2012l ) for details regarding justification of the mass 
range taken from each run. 

apply the DCsfh at 0.10 M© intervals using 

SFRDoU^, Ms) = SFRDint{z, Ms) x DC{z, Ms), (7) 

where SFRDohs{z, Ms) is what would be observable at 
a particular redshift, and SFRDint{z, Ms) is the intrinsic 
value in our simulation at the same redshift. Although we 
didn't indicate explicitly in the equation above, note that 
DC is also dependent on the values of limiting magnitude 
of the survey. 

This implementation results in decreasing the summed 
total SFRD (s olid magenta line w ith crosses) found in our 
previous work f Jaack s et aLll2012l ) by 71%, 73% and 77% at 
z — 6,7, 8, respectively, as shown by the solid green line. 
The most re cent observational da ta is shown by the black 
squares from lBouwens et aP (|201l[ ). Our results suggest that 
the DC of galaxies could reduce the observable fraction of 
SFRD down to only 20 — 30%, and it is necessary to correct 
for its effect in order to estimate the observable fraction of 
total SFRD in the Universe. 



5 CONCLUSIONS & DISCUSSIONS 
5.1 Main Conclusions 

Using cosmological SPH simulations, we examined the form 
of SFH of galaxies at z > 6, the duty cycle of star forma- 
tion and its impact on the current observations. Our major 
conclusions are as follows: 



• We find that, on average, the simulated SFH has an 
increasing form as a function of cosmic time, which can be 
characterized well by a power-law and an exponential func- 
tion (Fig. [2]). The best-fitting characteristic times are sum- 
marized in Tables [2] and [3] for different stellar mass range 
of galaxies and have characteristic time-scales of 70Myr to 
200 Myr for galaxies with stellar masses Ms ^ lO^M© to 
> 1O^°M0. We find that the power-law produces on aver- 
age slightly better fit to the mean SFH of galaxies, but the 
exponential form also works quite well. 

• The SFH of individual galaxies have bursts of star for- 
mation at z > 6, and these burst lead to a scatter in the 
relationship between Ms and Muv of AMuv ^ ±1 mag and 
a ^ 1 dex scatter in the Ms - SFR relation (Fig. g]). The 
measured val ues of duty cycl e refle ct the scatter in these 
relationships. iFeldmann et al. also supports the idea 
that star formation is inherently a stochastic process as they 
have implemented a model treating SF as a Poisson pro- 
cess within GMCs. iForero- Romero k, Diikstral (|2012h have 
shown that SF stochasticity can broaden the distribution 
of Lya equivalent width, and impact the observed scat- 
ter in the Muv — SFR relationship of Lya emitters with 
SFR < lO"^M0yr"^ These work support the idea that 
duty cycle can be a useful quantity to summarize the scat- 
ter in th e Ms — Muv — SFR re lationships. 

While iFinlator et al ] (|201lh find a very similar mean 
Ms — SFR relationship in their simulations they find a 
much smaller scatter which is on the order of 0.2 dex us- 
ing comparable resolution runs. They interpret this, along 
with their smoothly rising SFH, as evidence for smooth ac- 
cretion rather than merger-driven growth. In contrast, our 
^ 1 dex scatter suggests that we can not rule out the domi- 
nance of merger-driven growth for high-z galaxies, especially 
at z > 6 when the merger rates are very high and gravita- 
tional growth of structure is driving the rapid increase in 
cosmic SFRD. 

• We measured the duty cycle of simulated galaxies us- 
ing two different methods: 1) DCsfh, based on the galaxy's 
star formation history; 2) DCm^^, based on the instan- 
taneous UV magnitude of the galaxy. Both methods give 
DCs at z = 6, 7, 8 that make a steep transition from zero 
to unity at Ms = lO'^ - 10^ M©. The effect of dust ex- 
tinction comes in when we convert the magnitude limit 
into threshold SFR (SFRth) using the relation between 
the two quantities in our simulations. The range in which 
the transition occurs is highly dependent on the amount 
of dust extinction assumed to be present. With the dust 
effect of Eb-v = 0.10, the transition mass-scale is about 
log(Ms/M0) = 8.1 - 8.4, but without the dust effect, the 
transition occurs at log(Ms/M0) = 7.6 - 7.9 (see Table |4] 
andO. This may be better characterized as a shift to a lower 
DC for a given Ms as fewer galaxies will be able to exceed 
either the SFH or Muv thresholds as Eb-v increases. With 
a greater value of Eb-v, the galaxy has to have a higher 
SFR (and hence higher Ms) in order to exceed the observ- 
able threshold. 

• The application of DCsfh including Eb-v = 0.10 to 
the GSMF results in lowering the faint-end of the GSMF 
by on average Aa ~ 0.26 at z = 6, 7, 8 (Fig. I10|) . which is 
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simila r to the results found at lower redshifts bv iLee et alJ 

dmil). 

• The application of the duty cycle to the simulated 
SFRD results in lowering the total SFRD by 71, 72, 77% at 
z — 6,7,8, respectively fFig. [TT]). This indicates that the 
current observations (even the deepest HST galaxy imaging 
surveys probing down to Muv — 18 mag) could be missing 
a significant fraction (> 70%) of star formation in the early 
Universe. 



5.2 H2 regulated star formation 

The details of the SFHs of simulated galaxies could be de- 
pendent on the SF model implemented in our simulations. 
Therefore the di screpancies tha t we fi nd between our re- 
sults and that of iFinlator et all (|201lh may simply be due 
to the differences between the adopt ed SF model. Cur- 
rently we use the " Pressure mo del" ( Sc have k, Dalla Vecchial 
I2OO8I : IChoi Nagaminell201Ql ) , in which stars are allowed to 
form once the gas densit y exceeds the thr eshold density of 
n^if = 0.6 cm~^ (see Naga mine et al.ll20ld for the justifica- 
tion of this value). Star particles are then generated based 
on the gas pressure of st ar-forming ^ as acc ording to the SF 
law matched to the local iKennicuttI (| 19981 ) relationship. 

There is observational evidence to suggest that star 
form ation in the Un i verse i s controlled by the fraction of 
H2 (iKennicutt et al.1 l2007l : iLerov et al.1 120081 : iBigiel et al.1 
l2008i r This could lead to a reduction in star formation at 
high redshift due to lower metallicity (|Krumholz et al ] l2009l : 
iKuhlen et al]l2012l ). and would alter the SFHs in our simu- 
lation. Therefore this work would need to be revisited with 
a new SF model based on H2 mass. We are currently in the 
process of implementi ng a new SF model ba sed on H2 mass 
based on the work by (|Krumholz et al.ll2009 ) in our simula- 
tions, and we will report the results in a subsequent paper 
(Thompson, et al., in preparation). 
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